closefile (openfile ("physics-tb-problem.fee", "w"))
open("physics-tb-problem.fee")

internal_electrode_radius = 1.5e-3
external_electrode_radius = 2.0e-3
offset = 0
potential_difference = 3000
magnetic_induction = 3
angle = 0
delta_t = (1e-16)-- seconds T ~ 6.666666666e-12

showconsole()

print("internal_electrode_radius = ", internal_electrode_radius)
print("external_electrode_radius = ", external_electrode_radius)
print("offset  = ", offset)
print("potential_difference = ", potential_difference)
print("magnetic_induction = ", magnetic_induction)

-- Setup internal electrode
ei_addboundprop("Internal electrode", -potential_difference, 0, 0, 0, 0)

ei_addnode (-internal_electrode_radius, 0)
ei_addnode (internal_electrode_radius, 0)
ei_addarc (-internal_electrode_radius, 0,
				   internal_electrode_radius,  0,
					 180, 1)
ei_addarc (internal_electrode_radius,  0,
           -internal_electrode_radius, 0,
           180, 1)

ei_selectarcsegment (0, internal_electrode_radius)
ei_setarcsegmentprop (1, "Internal electrode", 0, 0, "<None>")
ei_clearselected ()

ei_selectarcsegment (0, -internal_electrode_radius)
ei_setarcsegmentprop (1, "Internal electrode", 0, 0, "<None>")
ei_clearselected ()


-- Setup external electrode
ei_addboundprop("External electrode", 0, 0, 0, 0, 0)

ei_addnode (-external_electrode_radius - offset, 0)
ei_addnode (external_electrode_radius - offset, 0)
ei_addarc (-external_electrode_radius - offset, 0,
           external_electrode_radius - offset,  0,
           180, 1)
ei_addarc (external_electrode_radius - offset,  0,
					 -external_electrode_radius - offset, 0,
           180, 1)

ei_selectarcsegment (0-offset, external_electrode_radius)
ei_setarcsegmentprop (1, "External electrode", 0, 0, "<None>")
ei_clearselected ()

ei_selectarcsegment (0-offset, -external_electrode_radius)
ei_setarcsegmentprop (1, "External electrode", 0, 0, "<None>")
ei_clearselected ()

-- "Space" inside internal electrode.
ei_addblocklabel (0,0)
ei_selectlabel (0,0)
ei_setblockprop("<No Mesh>", 1, 1, 0)
ei_clearselected()


-- Space between electrodes.
ei_getmaterial("Air")

ei_addblocklabel (-(internal_electrode_radius + (external_electrode_radius - internal_electrode_radius)/2 - offset), 0)
ei_selectlabel (-(internal_electrode_radius + (external_electrode_radius - internal_electrode_radius)/2 - offset), 0)
ei_setblockprop("Air", 0, 1e-5, 0)


-- Zoom to see everything
ei_zoomnatural()


-- Prepare problem solution
ei_probdef ("meters", "planar", 1.e-16)
ei_createmesh()

-- Solve problem
ei_analyze(0)
ei_loadsolution()

-- Postrocess problem
--showconsole()
clearconsole()

-- Electron parameters
m_e = 9.1093829140e-31   -- electron mass
q_e = -1.60217656530e-19 -- electron charge
n = q_e/m_e
w = n * magnetic_induction

--Initial conditions
--x0 = (internal_electrode_radius+1.e-8)*cos(angle)
--y0 = (internal_electrode_radius+1.e-8)*sin(angle)
x0 = (internal_electrode_radius)*cos(angle)
y0 = (internal_electrode_radius)*sin(angle)
dx0 = 0
dy0 = 0
Ex = 0
Ey = 0

-- Electron posotion


function x(t)
	 return (-(sin(t*w)*(n*Ey-dx0*w))/w^2+(t*n*Ey)/w+(n*Ex+w^2*x0+dy0*w)/w^2-(cos(t*w)*(n*Ex+dy0*w))/w^2)
end

function y(t)
	 return ((n*Ey+w^2*y0-dx0*w)/w^2-(cos(t*w)*(n*Ey-dx0*w))/w^2+(sin(t*w)*(n*Ex+dy0*w))/w^2-(t*n*Ex)/w)
end

function update_E (t)
	 if (t ~= 0) then
			_,_,_,Ex,Ey,_,_,_ = eo_getpointvalues (x(t), y(t))
   else
      _,_,_,Ex,Ey,_,_,_ = eo_getpointvalues (x0, y0)
	 end
	 
	 --print_current_parameters (t)
	 --print (Ex)
	 --print (Ey)
	 
   if (Ex ~= nil) and (Ey ~= nil) then
      Ex = Ex
      Ey = Ey
			--Ex = Dx/ex
			--Ey = Dy/ey
			--print ("here")
			return 1
	 else
			--print ("error at ", t)
			return 0
	 end
end

function print_current_parameters (t)
	 print ("t = ", t)
	 print ("\tx = ", x(t))
   print ("\ty = ", y(t))
   print ("\tEx = ", Ex)
   print ("\tEy = ", Ey)
	 print ("=========")
end

t = 0
steps = 0
update_E(0)
xs = {}
ys = {}
--eo_addcontour (x0, y0)
for i=0,1e6,1 do
   t = i*delta_t
   eo_addcontour (x(t), y(t))
	 if ((x(t)^2 +y(t)^2) >= external_electrode_radius^2)
	 then
			break
	 else
			xs [i] = x(t)
			ys [i] = y(t)
			steps = steps + 1
			--print_current_parameters (t)
	 end
	 
	 if (update_E(t) ~= 1) then
			break
	 end
end

--eo_addcontour (xs[steps - 1], ys[steps - 1])
--print_current_parameters (t)
print (xs[steps-2])
print (ys[steps-2])
print ("=======")
print (xs[steps-1])
print (ys[steps-1])

handle = openfile ("out.dat", "w")

for i= 0,steps - 1,1 do
   write (handle, xs[i])
   write (handle, ';')
   write (handle, ys[i])
   write (handle, '\n')
end

closefile(handle)

--quit()
